format short e
echo on
% La fattorizzazione LU senza pivot non puo' essere eseguita su alcune matrici
% eseguiamola sulla matrice
Q=[0 1 0;0 0 1;1 0 0 ]
pause
[L,U]=gauss(Q)
pause
A=[1 2 5;-1 -2 -7;1 1 3]
pause
[L,U]=gauss(A)
pause
[L,U,P]=lu(Q)
pause
[L,U,P]=lu(A)
pause
% presenta problemi di instabilita'  numerica
A=[ 1e-17 1; 1 1]
pause
[L,U]=gauss(A)
% Il prodotto LU non e' A  infatti:
norm(L*U  - A)
pause
% inoltre la matrice A e' ben condizionata
cond(A)
% mentre L e U non lo sono
cond(U)
cond(L)
pause
% eseguiamo la fattorizzazione con pivot
[L,U,P]=lu(A)
pause
% controlliamo i fattori di condizionamento
cond(L)
cond(U)
cond(A)
pause
%% Consideriamo il sistema lineare
A=[1.002 1;1 0.998]
pause
b=[2.002;1.998]
pause
%% La soluzione di tale sistema e' data da:
%% l'istruzione \ risolve il sistema utilizzando la 
%% fattorizzazione lu con pivot
x=A\b
pause
%% errore rispetto al valore esatto
norm(x -ones(2,1),'inf')
pause
%% calcoliamo il residuo
r = A*x-b
pause
%% residuo relativo
norm(r,'inf')/norm('inf')
pause
% numero di condizione
cond(A,inf)
pause
%% Perturbiamo b 
 bt = b;
 bt(1) = bt(1)+0.0001;
pause
%% errore relativo sui dati di input
norm(b-bt,'inf')/norm(b,'inf')
pause
%% ci aspettiamo che la soluzione del sistema perturbato
%% disti dalla soluzione del sistema di partenza di:
cond(A,'inf')*norm(b-bt,'inf')/norm(b,'inf')
pause
%% La soluzione del sistema perturbato e':
 xt=A\bt
 pause
% Ora possiamo confrontare le due soluzioni e valutare
%l'errore relativo:
[x xt]
pause
norm(x-xt)/norm(x)
pause
help gallery
pause
[A,b]=gallery('wilk',3)
pause
cond(A)
pause
x = A\b
pause
% residuo relativo
norm(A*x-b,'inf')/norm(b,'inf')
pause
% maggiorazione dell'errore
cond(A,'inf')*norm(A*x-b,'inf')/norm(b1,'inf')


%
%Ci sono sistemi che non sono mal condizionati ma per i
%quali il numero di condizione della matrice e' grande.
%Il condizionamento di una matrice pu\`o essere modificato dallo SCALING.
%Con tale metodo si va alla ricerca di due matrici diagonali D1 e D2
% in modo da minimizzare K(D1AD2).
pause
A
pause
D1=diag([1 1 1e10])
pause
D1*A
pause
cond(A,'inf')
pause
cond(D1*A,'inf')
pause
D2=diag([1e10 1 1])
pause
A*D2
pause
cond(A*D2,'inf')
pause
cond(D1*A*D2,'inf')
pause
[A,b]=gallery('wilk',4)
pause
cond(A)
pause
x=A\b
pause
% maggiorazione dell'errore
cond(A,'inf')*norm(A*x-b,'inf')/norm(b,'inf')
% soluzione teorica nota tutti 1
xt = ones(4,1)
pause
% termine noto
b1=A*xt
pause
% risolviamo il sistema
x1=A\b1
pause
% errore relativo
norm(x1-xt,'inf')
pause
A
pause
% matrici di scaling
D1 = diag([1e5 1 1 1])
pause
D2 = diag([1 1 1 1e5])
pause
D1*A
pause
D1*A*D2
pause
cond(D1*A*D2,1)
pause
% esempio con matrice random
A=rand(100);
pause
cond(A)
pause
[L,U]=gauss(A);
pause
xt = rand(100,1);
b=A*xt;
pause
type soll
pause
type solu
pause
y=soll(L,b);
pause
x=solu(U,y);
pause
norm(x-xt,'inf')/norm(xt,'inf')
pause
cond(L)
pause
cond(U)
pause
cond(A)
pause
% utilizziamo il pivot
[L,U,P]=lu(A);
pause
y=soll(L,P*b);
pause
xp=solu(U,y);
pause
norm(xp-xt,'inf')/norm(xt,'inf')
pause
cond(L)
pause
cond(U)
pause
cond(A)
pause
% altra matrice
A = gallery('prolate',20);
pause
cond(A)
pause
[L,U]=gauss(A);
pause
xt = rand(20,1);
b=A*xt;
pause
y=soll(L,b);
pause
x=solu(U,y);
pause
norm(x-xt,'inf')/norm(xt)
pause
cond(L)
pause
cond(U)
pause
cond(A)
pause
% utilizziamo il pivot
[L,U,P]=lu(A);
pause
y=soll(L,P*b);
pause
xp=solu(U,y);
pause
norm(xp-xt,'inf')/norm(xt,'inf')
pause
cond(L)
pause
cond(U)
pause
cond(A)
pause